In this dataset, we have FEV of 606 children, between the ages of 6 and 17, were measured. The dataset also provides additional information on these children: their age, their height, their gender and, most importantly, whether the child is a smoker or a non-smoker.
age: This feature is the age of the person in numeric digits. height: It is the height of the patient in inches. gender: Male or Female represented as ‘m’ and ‘f’ respectively smoking: Whether person/child smokes or not in 0 and 1’s. fev: The fev, which is an acronym for forced expiratory volume, is a measure of how much air a person can exhale (in litres) during a forced breath.This is the actual feature we will create model for.
Before loading the dataset I will load all the libraries.
require(tidyverse)
require(magrittr)
require(R2jags)
require(mcmcse)
require(bayesplot)
require(TeachingDemos)
Reading the dataset in FEVdata object.
FEVdata <- read.table(file="fev.txt",header=T, sep="\t")
head(FEVdata)
## age fev height gender smoking
## 1 9 1.708 57.0 f 0
## 2 8 1.724 67.5 f 0
## 3 7 1.720 54.5 f 0
## 4 9 1.558 53.0 m 0
## 5 9 1.895 57.0 m 0
## 6 8 2.336 61.0 f 0
Before doing anything lets just convert gender to 0 or 1’s. 0 repersents male and 1 represent female
FEVdata$gender = ifelse(FEVdata$gender=="m", 1, 0) # keep in mind we are replacing the original gender from character to 0 and 1s
Plotting between age and fev:
plot(x=FEVdata$age,y=FEVdata$fev,xlab = "Age",ylab = "FEV", col= c("darkgreen"))
We realized that fev increases linearly with the age.
Plotting between height and fev:
plot(x=FEVdata$height,y=FEVdata$fev,xlab = "Height",ylab = "FEV", col= c("darkblue"))
We realized that fev increases linearly with the height.
Plotting between gender and fev:
plot(x=FEVdata$gender,y=FEVdata$fev,xlab = "Female - Male",ylab = "FEV", col= c("darkgreen"))
We realized that male generally have more fev than female.
Plotting between smoking and fev
plot(x=FEVdata$smoking,y=FEVdata$fev,xlab = " Non-Smoking - Smoking",ylab = "FEV", col= c("darkblue"))
From above plot we didn’t get anything for now. Lets see we can conclude something here or not.
si = 0
ni = 0
smokers = 0
non_smokers = 0
n = length(FEVdata$smoking)
for(i in 1:n) {
if(FEVdata$smoking[i] == 1)
{
smokers[si] = FEVdata$fev[i]
si = si + 1
}
else
{
non_smokers[ni] = FEVdata$fev[i]
ni = ni + 1
}
}
mean_of_smoker = mean( smokers)
mean_of_nonsmoker = mean(non_smokers)
cat("Mean of Smoker:",mean_of_smoker,fill = TRUE)
## Mean of Smoker: 3.261483
cat("Mean of non-Smoker:", mean_of_nonsmoker)
## Mean of non-Smoker: 2.636331
We can see after calculating the mean that mean of smoker is higher than mean of non-smoker. It means that smoker has higher tendency of having higher fev
Now we plot the distribution of fev
hist(FEVdata$fev,breaks = 50,xlab = "FEV Distribution")
we can see from above that most children’s fev stays between 2-3.
We observe that fev follows normal distribution so we try to develop its distribution.
We will create a model in jags and analyze it. Lets only include Age and Smoke features and analyze them.
\[y \sim N(\mu,\tau) \\ \mu= \beta_1 + \beta_2(Age)+\beta_3(Smoke) \\ Prior\ Distribution \ of \ parameters \\ \beta_1 \sim N(0,0.001) \\ \beta_2 \sim N(0,0.001) \\ \beta_3 \sim N(0,0.001) \\ \tau \sim G(0.001,0.001)\]
model1.jags <- function() {
# Likelihood
for(i in 1:n){
fev[i] ~ dnorm(mu[i],tau)
mu[i] <- beta[1] + beta[2]*age[i] + beta[3]*smoke[i]
}
beta[1] ~ dnorm(0,0.001) # Diffuse prior
beta[2] ~ dnorm(0,0.001)
beta[3] ~ dnorm(0,0.001)
tau ~ dgamma(0.001,0.001) #gamma here is important
}
The Parameter \(\tau\) represents the variation in the data of fev. The Parameter \(\beta_1\) is linear transformation parameter used to balance out the mean of data. The Parameter \(\beta_2\) gives the weight to age variable, the larger this parameter means fev depends more on value to age. The Parameter \(\beta_3\) gives the weight to Smoke variable, the larger this parameter means FEV depends more on value to smoke.
From above we have the Distribution of the FEV and prior distribution we can run MCMC on the data and evaluate the parameters For initial point in Markov chain we set All \(\beta=0\) and \(\tau=1\) and we run 9000 thousands MCMC itrations with 3 chains
Lets prepare data for jags
# Preparing data for JAGS
n <- length(FEVdata$age)
fev <- FEVdata$fev
smoke <- FEVdata$smoking
age <- FEVdata$age
#Lets only include Age and Smoke features and analyze them
dat.jags <- list("n","fev","smoke", "age")
# Defining parameters of interest
mod.params <- c("beta","tau")
# Starting values
mod.inits <- function(){
list("tau" = 1, "beta" = c(0,0,0))
}
# Run JAGS
set.seed(1873829) # adding my matricula
mod1.fit <- jags(data = dat.jags, # DATA
model.file = model1.jags, inits = mod.inits, # MODEL
parameters.to.save = mod.params,
n.chains = 3, n.iter = 9000, n.burnin = 1000, n.thin=10)# MCMC
## module glm loaded
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 606
## Unobserved stochastic nodes: 4
## Total graph size: 1860
##
## Initializing model
mod1.fit
## Inference for Bugs model at "/var/folders/yy/6ydy4yyd44zclz1dgg2076lw0000gn/T//RtmphrLy0P/model639d37ea12ba.txt", fit using jags,
## 3 chains, each with 9000 iterations (first 1000 discarded), n.thin = 10
## n.sims = 2400 iterations saved
## mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat
## beta[1] 0.208 0.103 0.004 0.137 0.210 0.276 0.403 1.002
## beta[2] 0.248 0.010 0.228 0.241 0.247 0.255 0.267 1.002
## beta[3] -0.236 0.084 -0.402 -0.293 -0.236 -0.179 -0.068 1.002
## tau 3.135 0.182 2.788 3.010 3.133 3.256 3.497 1.002
## deviance 1030.038 2.928 1026.479 1027.888 1029.414 1031.413 1037.537 1.003
## n.eff
## beta[1] 1800
## beta[2] 1900
## beta[3] 1100
## tau 1100
## deviance 1100
##
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
##
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 4.3 and DIC = 1034.3
## DIC is an estimate of expected predictive error (lower deviance is better).
chainArray <- mod1.fit$BUGSoutput$sims.array
From above we get the beta-1 as 0.208, beta-2 value as 0.248, beta_3 as -0.236 and tau value as 1034.3. pD is 4.3 panalization value and DIC measure of goodness of fit throughtout different chains which has 1034.3.
lets get the estimates bayesian point and interval estimation of our model.
chainMat <- mod1.fit$BUGSoutput$sims.matrix
# Intervals
cred <- 0.95
cat("Here are pont estimates or means of different parameters:", fill = TRUE)
## Here are pont estimates or means of different parameters:
# Point estimates
(par.hat.jags <- colMeans(chainMat))
## beta[1] beta[2] beta[3] deviance tau
## 0.2076090 0.2475994 -0.2357305 1030.0375901 3.1349292
cat("", fill = TRUE)
cat("Here are equal tail intervals of our parameters we estimated:", fill = TRUE)
## Here are equal tail intervals of our parameters we estimated:
# Intervals
(par.ET.jags <- apply(chainMat, 2, quantile,
prob=c((1-cred)/2, 1-(1-cred)/2)))
## beta[1] beta[2] beta[3] deviance tau
## 2.5% 0.004133159 0.2283579 -0.40248049 1026.479 2.787757
## 97.5% 0.402949602 0.2673897 -0.06804479 1037.537 3.497325
cat("", fill = TRUE)
cat("Here are HPD intervals using coda of our model:", fill = TRUE)
## Here are HPD intervals using coda of our model:
# What about the HPD?
(par.HPD.jags <- coda::HPDinterval(as.mcmc(chainMat)))
## lower upper
## beta[1] 0.01211304 0.40702379
## beta[2] 0.22886315 0.26778785
## beta[3] -0.41606811 -0.08851262
## deviance 1026.11503822 1035.90043542
## tau 2.76242954 3.47113071
## attr(,"Probability")
## [1] 0.95
Lets do some dignostics using BayesPlot.
# Plots with BayesPlot
chainArray <- mod1.fit$BUGSoutput$sims.array
bayesplot::mcmc_combo(chainArray) # combo of density plot
We can see from bayesplot above that chains values of each parameters. It seems like our model is converged because different chains are exploring the same kind of area.
bayesplot::mcmc_acf(chainArray)
From above we see that our model is started from 1 and went quickly towards 0 or near zero for each parameters in chain.
Lets use coda now because it will provide us with geweke, gelmen and heidel diagnostics
coda.fit <- coda::as.mcmc(mod1.fit)
coda::geweke.diag(coda.fit) #
## [[1]]
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## beta[1] beta[2] beta[3] deviance tau
## -1.0803 1.0386 0.7099 0.1273 1.9995
##
##
## [[2]]
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## beta[1] beta[2] beta[3] deviance tau
## 0.4579 -0.1715 -1.2364 1.4817 1.5684
##
##
## [[3]]
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## beta[1] beta[2] beta[3] deviance tau
## 1.2742 -1.1209 0.3636 0.1248 1.2609
Geweke basically comparing mean of 1st 10% of chain to mean of last 50% of the chain if they are both kind of equal we surely converged. Lets plot it here:
coda::geweke.plot(coda.fit)
coda::gelman.diag(coda.fit)
## Potential scale reduction factors:
##
## Point est. Upper C.I.
## beta[1] 1 1.01
## beta[2] 1 1.01
## beta[3] 1 1.01
## deviance 1 1.00
## tau 1 1.00
##
## Multivariate psrf
##
## 1
coda::gelman.plot(coda.fit)
This is the potential scale reduction factor that I want to be near 1 or maybe under is better at the same time. The plots suggest to have the median (black line) under the 97.5% of quantile (red line).
lets now run the heidel diagnostic to anaylze that.
coda::heidel.diag(coda.fit)
## [[1]]
##
## Stationarity start p-value
## test iteration
## beta[1] passed 1 0.749
## beta[2] passed 1 0.753
## beta[3] passed 1 0.906
## deviance passed 1 0.699
## tau passed 81 0.128
##
## Halfwidth Mean Halfwidth
## test
## beta[1] passed 0.209 0.007196
## beta[2] passed 0.248 0.000706
## beta[3] passed -0.232 0.005729
## deviance passed 1030.057 0.206820
## tau passed 3.140 0.013204
##
## [[2]]
##
## Stationarity start p-value
## test iteration
## beta[1] passed 1 0.6800
## beta[2] passed 1 0.6779
## beta[3] passed 1 0.7395
## deviance passed 1 0.0801
## tau passed 1 0.4107
##
## Halfwidth Mean Halfwidth
## test
## beta[1] passed 0.211 0.007412
## beta[2] passed 0.247 0.000732
## beta[3] passed -0.234 0.005932
## deviance passed 1030.183 0.209045
## tau passed 3.125 0.012588
##
## [[3]]
##
## Stationarity start p-value
## test iteration
## beta[1] passed 1 0.646
## beta[2] passed 1 0.570
## beta[3] passed 1 0.825
## deviance passed 1 0.435
## tau passed 1 0.103
##
## Halfwidth Mean Halfwidth
## test
## beta[1] passed 0.203 0.007467
## beta[2] passed 0.248 0.000661
## beta[3] passed -0.241 0.005802
## deviance passed 1029.874 0.192069
## tau passed 3.137 0.012560
Heidel check the model in two steps first it check whether chain is stationary. From above we can see that chain seems to passed all stationary test. The second part basically checks that whether accuracy is good enough comparing it to the half width of chain.
lets create another model with slight modification and compare its DIC to one we have already created.
\[y \sim N(\mu,\tau) \\ \mu= \beta_1 + \beta_2(Age)+\beta_3(Smoke)+\beta_4(Age)(Smoke) \\ Prior\ Distribution \ of \ parameters \\ \beta_1 \sim N(0,0.001) \\ \beta_2 \sim N(0,0.001) \\ \beta_3 \sim N(0,0.001) \\ \beta_4 \sim N(0,0.001) \\ \tau \sim G(0.001,0.001)\]
In this model, the Parameter \(\beta_4\) gives the weight to product of smoke and age variables, the larger this parameters means fev depends more on product of Smoke and AgeHere is model in jags:
model2.jags <- function() {
# Likelihood
for(i in 1:n){
fev[i] ~ dnorm(mu[i],tau)
mu[i] <- beta[1] + beta[2]*age[i] + beta[3]*smoke[i] + beta[4]*age[i]*smoke[i]
}
beta[1] ~ dnorm(0,0.001) # Diffuse prior
beta[2] ~ dnorm(0,0.001)
beta[3] ~ dnorm(0,0.001)
beta[4] ~ dnorm(0,0.001)
tau ~ dgamma(0.001,0.001)
}
Lets run the model with same data we have already prepared in model 1.
# here we have to initialize the model beta_4 value equal to zero as well.
mod2.inits <- function(){
list("tau" = 1, "beta" = c(0,0,0,0))
}
set.seed(1873829)
mod2.fit <- jags(data = dat.jags, # DATA
model.file = model2.jags, inits = mod2.inits, # MODEL
parameters.to.save = mod.params,
n.chains = 3, n.iter = 9000, n.burnin = 1000, n.thin=10)# MCMC
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 606
## Unobserved stochastic nodes: 5
## Total graph size: 1882
##
## Initializing model
mod2.fit
## Inference for Bugs model at "/var/folders/yy/6ydy4yyd44zclz1dgg2076lw0000gn/T//RtmphrLy0P/model639d231dcd24.txt", fit using jags,
## 3 chains, each with 9000 iterations (first 1000 discarded), n.thin = 10
## n.sims = 2400 iterations saved
## mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat
## beta[1] 0.064 0.102 -0.130 -0.006 0.063 0.134 0.265 1.000
## beta[2] 0.262 0.010 0.242 0.255 0.262 0.269 0.281 1.001
## beta[3] 2.259 0.484 1.306 1.943 2.244 2.576 3.229 1.001
## beta[4] -0.193 0.037 -0.266 -0.216 -0.192 -0.168 -0.121 1.001
## tau 3.262 0.186 2.908 3.135 3.261 3.387 3.639 1.001
## deviance 1004.046 3.198 999.821 1001.664 1003.384 1005.712 1011.753 1.001
## n.eff
## beta[1] 2400
## beta[2] 2400
## beta[3] 2400
## beta[4] 2400
## tau 2400
## deviance 2400
##
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
##
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 5.1 and DIC = 1009.2
## DIC is an estimate of expected predictive error (lower deviance is better).
chainArray <- mod2.fit$BUGSoutput$sims.array
Through our DIC comparison we can see that DIC value of second model is slighly descreased it means that second model is better than the first model. I will not going to perform diagnostics in this model because intended parameters values will be just slighly changed.
Lets create a simulated data for our model estimates above to check how it is working so that we know that there is nothing wrong in the model itself as diagnostic check.
N = 606 # sample size
set.seed(1873829)
# Simulate the tries size
sim.smoke <- sample(c(0,1), replace=TRUE, size=N) # creating smoke
# Simulate the covariate (as you prefer)
sim.age <- sample(6:17, replace=TRUE, size=N)
beta = c(0.064,0.262 ,2.259,-0.193)
tau = 0.186 # standard diviation of tau
# Pick fixed values for the parameters of the model got from above model
sim.fev = 0
sim.mus = 0
for(i in 1:n){
sim.mus[i] <- beta[1] + beta[2]* sim.age[i] + beta[3]*sim.smoke[i]+ beta[4]*sim.age[i]*sim.smoke[i]
}
#Simulate response according to the model
sim.fev = rnorm(N, sim.mus, tau ) #predicting values based on our previous model.
sim.dat <- data.frame(age=sim.age, smoking= sim.smoke, fev=sim.fev)
head(sim.dat)
## age smoking fev
## 1 12 1 3.242332
## 2 6 0 1.124224
## 3 11 0 2.786389
## 4 11 1 3.149172
## 5 14 0 3.631531
## 6 11 0 3.057191
Same model there is a change of sim.model.jags.
sim.model.jags <- function() {
# Likelihood
for(i in 1:n){
fev[i] ~ dnorm(mu[i],tau)
mu[i] <- beta[1] + beta[2]* age[i] + beta[3]*smoke[i]+ beta[4]*age[i]*smoke[i]
}
beta[1] ~ dnorm(0, 0.001) # Diffuse prior
beta[2] ~ dnorm(0, 0.001)
beta[3] ~ dnorm(0, 0.001)
beta[4] ~ dnorm(0, 0.001)
tau ~ dgamma(0.001, 0.001)
}
Lets initialize the simulated model and run it.
# data that jags will use
sim.dat.jags <- list(n=N,fev=sim.fev,age=sim.age,smoke=sim.age)
# parameters of intrests
sim.mod.params <- c("beta","tau")
# Starting values
sim.mod.inits <- function(){
list("tau" = 1, "beta" = c(0,0,0,0))
}
set.seed(1873829)
sim.mod.fit <- jags(data = sim.dat.jags, # DATA
model.file = sim.model.jags, inits = sim.mod.inits, # MODEL
parameters.to.save = sim.mod.params,
n.chains = 3, n.iter = 9000, n.burnin = 1000, n.thin=10) # MCMC
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 606
## Unobserved stochastic nodes: 5
## Total graph size: 1874
##
## Initializing model
sim.mod.fit
## Inference for Bugs model at "/var/folders/yy/6ydy4yyd44zclz1dgg2076lw0000gn/T//RtmphrLy0P/model639d9606104.txt", fit using jags,
## 3 chains, each with 9000 iterations (first 1000 discarded), n.thin = 10
## n.sims = 2400 iterations saved
## mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
## beta[1] 1.439 0.189 1.079 1.309 1.437 1.570 1.811 1.000 2400
## beta[2] -0.397 22.969 -46.269 -15.937 -0.009 15.331 44.406 1.001 2400
## beta[3] 0.529 22.969 -44.259 -15.196 0.137 16.045 46.416 1.001 2400
## beta[4] 0.001 0.001 -0.002 0.000 0.001 0.002 0.004 1.000 2400
## tau 6.700 0.382 5.979 6.440 6.699 6.955 7.485 1.001 2400
## deviance 567.899 2.835 564.341 565.811 567.252 569.381 575.020 1.001 2400
##
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
##
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 4.0 and DIC = 571.9
## DIC is an estimate of expected predictive error (lower deviance is better).
It seems like from above that DIC value descreased quite a bit also model was able to get the interested parameter. It is performing exceptionally welll on simulated data that model requires per its implementation.
Just to be sure lets perform coda heidel test
coda.fit <- coda::as.mcmc(sim.mod.fit)
coda::heidel.diag(coda.fit)
## [[1]]
##
## Stationarity start p-value
## test iteration
## beta[1] passed 1 0.0662
## beta[2] passed 1 0.5688
## beta[3] passed 1 0.5677
## beta[4] passed 1 0.1500
## deviance passed 1 0.6474
## tau passed 1 0.7647
##
## Halfwidth Mean Halfwidth
## test
## beta[1] passed 1.44112 0.013128
## beta[2] failed -0.32697 1.732909
## beta[3] failed 0.45856 1.733017
## beta[4] passed 0.00114 0.000102
## deviance passed 567.92646 0.191424
## tau passed 6.70859 0.026958
##
## [[2]]
##
## Stationarity start p-value
## test iteration
## beta[1] passed 1 0.207
## beta[2] passed 1 0.700
## beta[3] passed 1 0.699
## beta[4] passed 1 0.377
## deviance passed 1 0.142
## tau passed 1 0.324
##
## Halfwidth Mean Halfwidth
## test
## beta[1] passed 1.44e+00 1.26e-02
## beta[2] failed 1.05e-01 1.61e+00
## beta[3] failed 2.70e-02 1.61e+00
## beta[4] passed 1.14e-03 8.57e-05
## deviance passed 5.68e+02 2.03e-01
## tau passed 6.69e+00 2.85e-02
##
## [[3]]
##
## Stationarity start p-value
## test iteration
## beta[1] passed 1 0.0549
## beta[2] passed 1 0.8709
## beta[3] passed 1 0.8692
## beta[4] passed 1 0.0518
## deviance passed 1 0.6644
## tau passed 1 0.2470
##
## Halfwidth Mean Halfwidth
## test
## beta[1] passed 1.43976 0.013474
## beta[2] failed -0.96874 1.575041
## beta[3] failed 1.10041 1.575051
## beta[4] passed 0.00114 0.000103
## deviance passed 567.84628 0.195409
## tau passed 6.69950 0.024182
Only halfwidth Mean test is failing. Lets to be sure we will perform one other diagnostic test.
coda::gelman.diag(coda.fit)
## Potential scale reduction factors:
##
## Point est. Upper C.I.
## beta[1] 1 1.01
## beta[2] 1 1.00
## beta[3] 1 1.00
## beta[4] 1 1.01
## deviance 1 1.00
## tau 1 1.00
##
## Multivariate psrf
##
## 1
coda::gelman.plot(coda.fit)
After performing the analysis we can say that our model was able to get the model parameters correctly based on simulated data.
We will create a linear model from the data and try to fit the parameters. Linear regressin is good choice for frequentistic approach. For this analysis we will get only two features age and smoke and do the same analysis we used in baysian.
smoke<-FEVdata$smoking
fev<- FEVdata$fev
age<- FEVdata$age
n <- length(fev)
lr <- lm(fev~age+smoke+age:smoke) # linear regression
# should display summary of the model
summary(lr)
##
## Call:
## lm(formula = fev ~ age + smoke + age:smoke)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.39397 -0.35731 -0.02474 0.31752 2.01363
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.06470 0.10087 0.641 0.522
## age 0.26214 0.01000 26.213 < 2e-16 ***
## smoke 2.26181 0.48397 4.673 3.66e-06 ***
## age:smoke -0.19292 0.03685 -5.236 2.28e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5536 on 602 degrees of freedom
## Multiple R-squared: 0.5569, Adjusted R-squared: 0.5547
## F-statistic: 252.2 on 3 and 602 DF, p-value: < 2.2e-16
beta <- coef(lr)
X <- cbind(rep(1,n),age,smoke,age*smoke)
# building getting data from model itself
rstudsmoke <- rstudent(lr)[smoke==1]
rstudnonsmoke <- rstudent(lr)[smoke==0]
We can observer that we get estimates of the parameters. Now we analyze our linear line with the dataset from which the model is generated. Fitted is a generic function which extracts fitted values from objects returned by modeling functions. Residuals extracts model residuals from objects returned by modeling functions.
par(mfrow=c(2,2))
plot(fitted(lr),resid(lr),xlab="Fitted Values", ylab="Residuals",main="", )
lines(lowess(fitted(lr),resid(lr)),col= c("red"))
qqnorm(resid(lr),main="")
qqline(resid(lr), col=c("red"))
Now we analyze our linear line with the dataset for smokers only
par(mfrow=c(2,2))
plot(fitted(lr)[smoke==1],rstudsmoke,xlab="Fitted Values", ylab="Studentized Residuals",main="Smokers \n Residual Plot")
lines(lowess(fitted(lr)[smoke==1],rstudsmoke), col= c("red"))
qqnorm(rstudsmoke,main="Smokers \n Normal Plot")
qqline(rstudsmoke,col= c("red"))
Now we analyize our linear line with the dataset for non-smokers only
par(mfrow=c(2,2))
plot(fitted(lr)[smoke==0],rstudnonsmoke,xlab="Fitted Values", ylab="Studentized Residuals",main="Nonsmokers \n Residual Plot")
lines(lowess(fitted(lr)[smoke==0],rstudnonsmoke), col= c("red"))
qqnorm(rstudnonsmoke,main="Nonsmokers \n Normal Plot")
qqline(rstudnonsmoke,col= c("red"))
Now we analyize our linear line with the dataset Age wise.
plot(age, rstudent(lr),xlab="Age", ylab="Residuals",main="")
lines(lowess(age,rstudent(lr)), col= c("red"))
Finally we are going to analyze our frequentistic model with model 2(More Accurate) which has less DIC value.
summary(lr)
##
## Call:
## lm(formula = fev ~ age + smoke + age:smoke)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.39397 -0.35731 -0.02474 0.31752 2.01363
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.06470 0.10087 0.641 0.522
## age 0.26214 0.01000 26.213 < 2e-16 ***
## smoke 2.26181 0.48397 4.673 3.66e-06 ***
## age:smoke -0.19292 0.03685 -5.236 2.28e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5536 on 602 degrees of freedom
## Multiple R-squared: 0.5569, Adjusted R-squared: 0.5547
## F-statistic: 252.2 on 3 and 602 DF, p-value: < 2.2e-16
mod2.fit
## Inference for Bugs model at "/var/folders/yy/6ydy4yyd44zclz1dgg2076lw0000gn/T//RtmphrLy0P/model639d231dcd24.txt", fit using jags,
## 3 chains, each with 9000 iterations (first 1000 discarded), n.thin = 10
## n.sims = 2400 iterations saved
## mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat
## beta[1] 0.064 0.102 -0.130 -0.006 0.063 0.134 0.265 1.000
## beta[2] 0.262 0.010 0.242 0.255 0.262 0.269 0.281 1.001
## beta[3] 2.259 0.484 1.306 1.943 2.244 2.576 3.229 1.001
## beta[4] -0.193 0.037 -0.266 -0.216 -0.192 -0.168 -0.121 1.001
## tau 3.262 0.186 2.908 3.135 3.261 3.387 3.639 1.001
## deviance 1004.046 3.198 999.821 1001.664 1003.384 1005.712 1011.753 1.001
## n.eff
## beta[1] 2400
## beta[2] 2400
## beta[3] 2400
## beta[4] 2400
## tau 2400
## deviance 2400
##
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
##
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 5.1 and DIC = 1009.2
## DIC is an estimate of expected predictive error (lower deviance is better).
We see that both models have given the output parameters and parameters in frequentist models are almost equal to mean of distribution of parameters in Baysean model. In baysean model we have fev as Normal distribution with variance as τ so baysian model also generates the variance of the fev in comparison to exact parameters in frequentist approach.